OHI British Columbia | OHI Science | Citation policy
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(sf)
library(raster)
dir_ohibc <- '~/github/ohibc'
dir_calc <- file.path(dir_ohibc, 'calc_ohibc')
dir_master <- file.path(dir_calc, 'master')
dir_figs <- file.path(dir_ohibc, 'ms_figures')
### provenance tracking
# library(provRmd); prov_setup()
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
scores_all <- read_csv(file.path(dir_calc, 'scores_all.csv'))
goals_all <- scores_all$goal %>% unique()ohibc_sf <- read_sf(file.path(dir_ohibc, 'prep/_spatial/ohibc_rgns_unclipped.shp'))
if(!file.exists(file.path(dir_figs, 'etopo_rast.tif'))) {
### Use ETOPO raster for bathy and topo shading
etopo_file <- file.path(dir_M, 'git-annex/bcprep/_raw_data',
'etopo-bathymetry/d2015/ETOPO2v2_bathymetry/ETOPO2v2c.tif')
etopo_raw <- raster(etopo_file)
crs(etopo_raw) <- '+init=epsg:4326'
extent(etopo_raw) <- c(-180, 180, -90, 90)
bbox <- st_bbox(ohibc_sf)[c(1, 3, 2, 4)] + c(-25000, +25000, -25000, +25000)
etopo <- etopo_raw %>%
crop(extent(-140, -110, 45, 65)) %>%
projectRaster(crs = '+init=epsg:3005', res = 1000) %>%
crop(extent(bbox))
bcland_sf <- read_sf(file.path(dir_ohibc, 'prep/_spatial/ohibc_land.shp')) %>%
st_set_crs(3005) %>%
st_crop(etopo)
bcland_rast <- etopo %>%
mask(bcland_sf %>% as('Spatial'))
ohibc_rast <- etopo %>%
mask(ohibc_sf %>% as('Spatial'))
values(etopo)[values(etopo) > -5 &
is.na(values(bcland_rast)) &
!is.na(values(ohibc_rast))] <- -5
### values close (or above) zero that aren't "land" but are in OHIBC,
### set to a nice depth
writeRaster(etopo, file.path(dir_figs, 'etopo_rast.tif'), overwrite = TRUE)
} else {
etopo <- raster(file.path(dir_figs, 'etopo_rast.tif'))
crs(etopo) <- '+init=epsg:3005'
}
slope_raster <- terrain(etopo, opt = 'slope')
aspect_raster <- terrain(etopo, opt = 'aspect')
hill_raster <- hillShade(slope_raster, aspect_raster, 40, 270)
hill_df <- rasterToPoints(hill_raster) %>%
as.data.frame() %>%
setNames(c('x', 'y', 'hill')) %>%
mutate(hill = ifelse(hill > .5, hill * 2 - 1, 0))
### just emphasize the steep hills
bcland_sf <- read_sf(file.path(dir_ohibc, 'prep/_spatial/ohibc_land.shp')) %>%
st_set_crs(3005) %>%
st_crop(etopo)
etopo_df <- etopo %>%
rasterToPoints() %>%
as.data.frame() %>%
setNames(c('x', 'y', 'elev'))
elev_r <- range(etopo_df$elev, na.rm = TRUE)
# basemap_cols <- c(rev(brewer.pal(5, 'Blues')), terrain.colors(5))
basemap_cols <- c(rev(brewer.pal(5, 'Blues')), 'forestgreen', 'green4', '#CCCCCC', '#EEEEEE', '#FFFFFF')
### set up values for where colors land, so zero matches up
val_0 <- (-elev_r[1]) / (elev_r[2] - elev_r[1])
basemap_vals <- c(seq(0, val_0, length.out = 5), seq(val_0 + .001, 1, length.out = 5))
basemap <- ggplot() +
ggtheme_map() +
geom_raster(data = etopo_df, aes(x, y, fill = elev), alpha = .5, show.legend = FALSE) +
geom_raster(data = hill_df, aes(x, y, alpha = hill), color = 'black', show.legend = FALSE) +
scale_fill_gradientn(colors = basemap_cols,
values = basemap_vals) +
scale_alpha_continuous(range = c(0, .5)) +
# geom_sf(data = bcland_sf, alpha = .5, color = 'grey60', fill = NA, size = .25) +
geom_sf(data = ohibc_sf, alpha = .5, color = 'red', fill = NA, size = .5) +
coord_sf(datum = NA) +
scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))scores <- scores_all %>%
filter(year == 2016) %>%
filter(dimension == 'score') %>%
complete(goal = goals_all, nesting(region_id))
### make sure every region has every goal, NA or otherwise
source('fig_fxns.R')
overall_flower <- plot_flower(scores %>% filter(region_id == 0))
# print(overall_flower)
rgn_flowers <- vector('list', length = 8)
rgn_names <- get_rgn_names()
for(rgn in 1:8) { ### rgn <- 1
rgn_name <- sprintf('%s (%s)',
rgn_names$rgn_name[rgn_names$rgn_id == rgn],
rgn_names$rgn_code[rgn_names$rgn_id == rgn])
rgn_name <- rgn_name %>% str_replace('est V|orth V', '. V')
### abbreviate West and North Vancouver Island (not North Coast)
rgn_name <- rgn_name %>% str_replace('Island', 'Is.')
### abbreviate Island for those regions
rgn_flowers[[rgn]] <- plot_flower(scores %>% filter(region_id == rgn),
score_size = 6,
plot_clean = TRUE,
rgn_name = rgn_name)
# print(rgn_flowers[[rgn]])
}library(cowplot)
sm_fl_size <- 0.21
lg_fl_size <- 0.60
draw_plot_ctr <- function(grob, x_ctr, y_ctr, size) {
x <- x_ctr - .5 * size; y <- y_ctr - .5 * size
draw_plot(grob, x, y, size, size)
}
x <- ggdraw() +
draw_plot(basemap, 0, 0, 1, 1) +
draw_plot_ctr(overall_flower, .30, .30, lg_fl_size) +
draw_plot_ctr(rgn_flowers[[1]], 0.55, 0.80, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[2]], 0.37, 0.87, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[3]], 0.65, 0.58, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[4]], 0.82, 0.54, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[5]], 0.70, 0.17, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[6]], 0.88, 0.35, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[7]], 0.20, 0.65, sm_fl_size) +
draw_plot_ctr(rgn_flowers[[8]], 0.72, 0.77, sm_fl_size)
ggsave('flower_plot_map.png', height = 9, width = 10, dpi = 300)